#ifndef BACKDIFFERFORMULAR_HPP
#define BACKDIFFERFORMULAR_HPP

#include <iostream>
#include <vector>
#include <math.h>
#include <fstream>
#include "norm.hpp"
#include "Function.hpp"

using namespace std;

class BackDifferFormula : public TimeIntegrator
{
private:
    int n_steps, p;
    double T, mu;
    vector<double> value;
    vector<vector<double>> result;
    Norm norm;
public:
    BackDifferFormula() {}
    BackDifferFormula(int n_, int _p, double time, vector<double> Val) : n_steps(n_), p(_p), T(time), value(Val) {}
    ~BackDifferFormula() {}
    void solver();
    void error();
};

void BackDifferFormula::solver()
{
    ofstream outfile;
    outfile.open("output.m");
    outfile << "points = [";
    double k = 1.0 * T / n_steps;
    vector<double> U = value;
    outfile << U << endl;
    result.push_back(U);
    for (int i = 0; i < p - 1; i++)
    {
        auto y1 = f(U, i * k);
        auto y2 = f(U + (0.5 * k) * y1, i * k + 0.5 * k);
        auto y3 = f(U + (0.5 * k) * y2, i * k + 0.5 * k);
        auto y4 = f(U + k * y3, (i + 1) * k);
        U = U + (1.0 * k / 6) * (y1 + 2.0 * y2 + 2.0 * y3 + y4);
        outfile << U << endl;
        result.push_back(U);
    }
    if (p == 1)
    {
        for (int i = 1; i < n_steps; i++)
        {
            auto tmp = U = U + k * f(U, k * i);
            U = result[i - 1] + k * f(tmp, i * k);
            outfile << U << endl;
            result.push_back(U);
        }
    }
    else if (p == 2)
    {
        for (int i = 1; i < n_steps; i++)
        {
            auto tmp = result[i] + k * (-0.5 * f(result[i - 1], (i - 1) * k) + 1.5 * f(result[i], i * k));
            U = 4.0 / 3 * result[i] - 1.0 / 3 * result[i - 1] + k * 2.0 / 3 * f(tmp, (i + 1) * k);
            outfile << U << endl;
            result.push_back(U);
        }
    }
    else if (p == 3)
    {
        for (int i = 2; i < n_steps; i++)
        {
            auto tmp = result[i] + k * (5.0 / 12 * f(result[i - 2], (i - 2) * k) - 4.0 / 3 * f(result[i - 1], (i - 1) * k) + 23.0 / 12 * f(result[i], i * k));
            U = 18.0 / 11 * result[i] - 9.0 / 11 * result[i - 1] + 2.0 / 11 * result[i - 2] + k * 6.0 / 11 * f(tmp, (i + 1) * k);
            outfile << U << endl;
            result.push_back(U);
        }
    }
    else if (p == 4)
    {
        for (int i = 3; i < n_steps; i++)
        {
            auto tmp = result[i] + k * (-9.0 / 24 * f(result[i - 3], (i - 3) * k) + 37.0 / 24 * f(result[i - 2], (i - 2) * k) - 59.0 / 24 * f(result[i - 1], (i - 1) * k) + 55.0 / 24 * f(result[i], i * k));
            U = 48.0 / 25 * result[i] - 36.0 / 25 * result[i - 1] + 16.0 / 25 * result[i - 2] - 3.0 / 25 * result[i - 3] + k * 12.0 / 25 * f(tmp, (i + 1) * k);
            outfile << U << endl;
            result.push_back(U);
        }
    }
    else
    {
        cout << "No such p for adamsmoulton." << endl;
    }
    outfile << "];" << endl;
    outfile << "scatter(points(:,1), points(:,2),\".\");" << endl;
    outfile.close();
}


void BackDifferFormula::error()
{
    // cout << result[n_steps] << endl;
    // cout << result[0] << endl;
    double error = norm.Inftynorm(result[n_steps] - result[0]);
    cout << error << endl;
}

TimeIntegrator* createBackDifferFormula(int n_, int p_, double T_, vector<double> value_) { return new BackDifferFormula(n_, p_, T_, value_); }
#endif
